home *** CD-ROM | disk | FTP | other *** search
/ Amiga Format CD 46 / Amiga Format CD46 (1999-10-20)(Future Publishing)(GB)[!][issue 1999-12].iso / -in_the_mag- / reader_requests / scilab / tests / dassldasrt.tst < prev    next >
Text File  |  1999-09-16  |  8KB  |  222 lines
  1. //DASSL
  2. // PROBLEM 1..   LINEAR DIFFERENTIAL/ALGEBRAIC SYSTEM
  3. //
  4. //X1DOT + 10.0*X1 = 0  
  5. //X1 + X2 = 1
  6. //X1(0) = 1.0, X2(0) = 0.0
  7. //
  8. t=1:10;t0=0;y0=[1;0];y0d=[-10;0];
  9. info=list([],0,[],[],[],0,0);
  10. //    Calling Scilab functions
  11. deff('[r,ires]=dres1(t,y,ydot)','r=[ydot(1)+10*y(1);y(2)+y(1)-1];ires=0')
  12. comp(dres1);
  13. deff('pd=djac1(t,y,ydot,cj)','pd=[cj+10,0;1,1]')
  14. comp(djac1);
  15. //   scilab function, without jacobian
  16. yy0=dassl([y0,y0d],t0,t,dres1,info);
  17. //   scilab functions, with jacobian
  18. yy1=dassl([y0,y0d],t0,t,dres1,djac1,info);
  19. // fortran routine, without jocabian
  20. yy2=dassl([y0,y0d],t0,t,'dres1',info);   //=yy0
  21. if norm(yy2-yy0,1)>1E-5 then pause,end
  22. // fortran routines, with jacobian
  23. yy3=dassl([y0,y0d],t0,t,'dres1','djac1',info);  //=yy1
  24. if norm(yy3-yy1,1)>1E-5 then pause,end
  25. yy3bis=dassl([y0,y0d],t0,t,'dres1',djac1,info); 
  26. // call fortran dres1 and scilab's djac1
  27. yy3ter=dassl([y0,y0d],t0,t,dres1,'djac1',info);
  28. //
  29. // with specific atol and rtol parameters
  30. atol=1.d-6;rtol=0;
  31. yy4=dassl([y0,y0d],t0,t,atol,rtol,dres1,info);
  32. yy5=dassl([y0,y0d],t0,t,atol,rtol,'dres1',info); //=yy4
  33. if norm(yy5-yy4,1)>1E-9 then pause,end
  34. yy6=dassl([y0,y0d],t0,t,atol,rtol,dres1,djac1,info); 
  35. yy7=dassl([y0,y0d],t0,t,atol,rtol,'dres1','djac1',info); //==yy6
  36. if norm(yy7-yy6,1)>1E-12 then pause,end
  37. //    
  38. //   Testing E xdot - A x=0
  39. //   x(0)=x0;   xdot(0)=xd0
  40. rand('seed',0);
  41. nx=5;
  42. E=rand(nx,1)*rand(1,nx);A=rand(nx,nx);
  43. //         Index 1
  44. [Si,Pi,Di,o]=penlaur(E,A);pp=Si*E;[q,m]=fullrf(pp);x0=q(:,1);x0d=pinv(E)*A*x0;
  45. deff('[r,ires]=g(t,x,xdot)','r=E*xdot-A*x;ires=0');comp(g);
  46. t=[1,2,3];t0=0;info=list([],0,[],[],[],0,0);
  47. x=dassl([x0,x0d],t0,t,g,info);x1=x(2:nx+1,:);
  48. if norm(pp*x1-x1,1)>1.d-5 then pause,end
  49. //x(4) goes through 1 at  t=1.5409711;
  50. t=1.5409711;ww=dassl([x0,x0d],t0,t,g,info);
  51. if abs(ww(5)-1)>0.001 then pause,end
  52. deff('[rt]=surface(t,y,yd)','rt=y(4)-1');nbsurf=1;
  53. [yyy,nnn]=dasrt([x0,x0d],t0,t,g,nbsurf,surface,info);
  54. deff('pd=j(t,y,ydot,cj)','pd=cj*E-A');comp(j);
  55. x=dassl([x0,x0d],t0,t,g,j,info);x2=x(2:nx+1,1);
  56. if norm(x2-ww(2:nx+1,1),1)>0.0001 then pause,end
  57. [yyy1,nnn]=dasrt([x0,x0d],t0,t,g,j,nbsurf,surface,info);
  58. //x0d is not known:
  59. x0d=ones(x0);info(7)=1;
  60. x=dassl([x0,x0d],t0,t,g,info);
  61. xn=dassl([x0,x0d],t0,t,g,j,info);
  62. if norm(x-xn,1)>0.00001 then pause,end
  63.  
  64.  
  65. //PROBLEM 2..
  66.  
  67. info=list([],0,[],[],[],0,0);
  68. y0=zeros(25,1);y0(1)=1;
  69. delta=0*y0;
  70. //link('dres2.o','dres2');
  71. //y0d=fort('dres2',0,1,'d',y0,2,'d',delta,3,'d',0,5,'i',0,6,'d',0,7,'d','out',[25,1],4,'d');
  72. y0d=zeros(y0);y0d(1)=-2;y0d(2)=1;y0d(6)=1;
  73. t0=0;t=[0.01,0.1,1,10,100];
  74. rtol=0;atol=1.d-6;
  75. y=dassl([y0,y0d],t0,t,atol,rtol,'dres2',info);
  76.  
  77. //                 DASRT
  78. // 
  79. //C-----------------------------------------------------------------------
  80. //C First problem.
  81. //C The initial value problem is..
  82. //C   DY/DT = ((2*LOG(Y) + 8)/T - 5)*Y,  Y(1) = 1,  1 .LE. T .LE. 6
  83. //C The solution is  Y(T) = EXP(-T**2 + 5*T - 4), YPRIME(1) = 3
  84. //C The two root functions are..
  85. //C   G1 = ((2*LOG(Y)+8)/T - 5)*Y (= DY/DT)  (with root at T = 2.5),
  86. //C   G2 = LOG(Y) - 2.2491  (with roots at T = 2.47 and 2.53)
  87. //C-----------------------------------------------------------------------
  88. y0=1;t=2:6;t0=1;y0d=3;
  89. info=list([],0,[],[],[],0,0);
  90. atol=1.d-6;rtol=0;ng=2;
  91. [yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,'res1',ng,'gr1',info);
  92. if abs(nn(1)-2.47)>0.001 then pause,end
  93. y0=yy(2,2);y0d=yy(3,2);t0=nn(1);t=[3,4,5,6];
  94. [yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,'res1',ng,'gr1',info);
  95. if abs(nn(1)-2.5)>0.001 then pause,end
  96. y0=yy(2,1);y0d=yy(3,1);t0=nn(1);t=[3,4,5,6];
  97. [yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,'res1',ng,'gr1',info);
  98. if abs(nn(1)-2.53)>0.001 then pause,end
  99.  
  100. deff('[delta,ires]=res1(t,y,ydot)','ires=0;delta=ydot-((2*log(y)+8)/t-5)*y')
  101. deff('[rts]=gr1(t,y,yd)','rts=[((2*log(y)+8)/t-5)*y;log(y)-2.2491]')
  102. comp(res1);comp(gr1);
  103.  
  104. y0=1;t=2:6;t0=1;y0d=3;
  105. info=list([],0,[],[],[],0,0);
  106. atol=1.d-6;rtol=0;ng=2;
  107. [yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,res1,ng,gr1,info);
  108. if abs(nn(1)-2.47)>0.001 then pause,end
  109. y0=yy(2,2);y0d=yy(3,2);t0=nn(1);t=[3,4,5,6];
  110. [yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,res1,ng,gr1,info);
  111. if abs(nn(1)-2.5)>0.001 then pause,end
  112. y0=yy(2,1);y0d=yy(3,1);t0=nn(1);t=[3,4,5,6];
  113. [yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,res1,ng,gr1,info);
  114. if abs(nn(1)-2.53)>0.001 then pause,end
  115.  
  116. //C
  117. //C-----------------------------------------------------------------------
  118. //C Second problem (Van Der Pol oscillator).
  119. //C The initial value problem is..
  120. //C   DY1/DT = Y2,  DY2/DT = 100*(1 - Y1**2)*Y2 - Y1,
  121. //C   Y1(0) = 2,  Y2(0) = 0,  0 .LE. T .LE. 200
  122. //C   Y1PRIME(0) = 0, Y2PRIME(0) = -2
  123. //C The root function is  G = Y1.
  124. //C An analytic solution is not known, but the zeros of Y1 are known
  125. //C to 15 figures for purposes of checking the accuracy.
  126. //C-----------------------------------------------------------------------
  127. rtol=[1.d-6;1.d-6];atol=[1.d-6;1.d-4];
  128. t0=0;y0=[2;0];y0d=[0;-2];t=[20:20:200];ng=1;
  129. info=list([],0,[],[],[],0,0);
  130. [yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,'res2','jac2',ng,'gr2',info);
  131. if abs(nn(1)-81.163512)>0.001 then pause,end
  132.  
  133. deff('[delta,ires]=res2(t,y,ydot)',...
  134. 'ires=0;y1=y(1),y2=y(2),delta=[ydot-[y2;100*(1-y1*y1)*y2-y1]]')
  135. [yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,res2,'jac2',ng,'gr2',info);
  136. deff('J=jac2(t,y,ydot,c)','y1=y(1);y2=y(2);J=[c,-1;200*y1*y2+1,c-100*(1-y1*y1)]')
  137. [yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,res2,jac2,ng,'gr2',info);
  138. deff('s=gr2(t,y,yd)','s=y(1)')
  139. [yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,res2,jac2,ng,gr2,info);
  140.  
  141. //           Hot Restart
  142.  
  143. [yy,nn,hot]=dasrt([y0,y0d],t0,t,atol,rtol,'res2','jac2',ng,'gr2',info);
  144. t01=nn(1);t=100:20:200;[pp,qq]=size(yy);y01=yy(2:3,qq);y0d1=yy(3:4,qq);
  145. [yy,nn,hot]=dasrt([y01,y0d1],t01,t,atol,rtol,'res2','jac2',ng,'gr2',info,hot);
  146. if abs(nn(1)-162.57763)>0.001 then pause,end
  147.  
  148. //Test of Dynamic link (Require f77!)
  149. //         1 making the routines
  150. res22=[...
  151. '      SUBROUTINE RES22(T,Y,YDOT,DELTA,IRES,RPAR,IPAR)';
  152. '      IMPLICIT DOUBLE PRECISION (A-H,O-Z)';
  153. '      INTEGER NEQ';
  154. '      DIMENSION Y(*), YDOT(*), DELTA(*)';
  155. '      NEQ=2';
  156. '      CALL F2(NEQ,T,Y,DELTA)';
  157. '      DO 10 I = 1,NEQ';
  158. '         DELTA(I) = YDOT(I) - DELTA(I)';
  159. ' 10   CONTINUE';
  160. '      RETURN';
  161. '      END';
  162. '      SUBROUTINE F2 (NEQ, T, Y, YDOT)';
  163. '      IMPLICIT DOUBLE PRECISION (A-H,O-Z)';
  164. '      INTEGER NEQ';
  165. '      DOUBLE PRECISION T, Y, YDOT';
  166. '      DIMENSION Y(*), YDOT(*)';
  167. '      YDOT(1) = Y(2)';
  168. '      YDOT(2) = 100.0D0*(1.0D0 - Y(1)*Y(1))*Y(2) - Y(1)';
  169. '      RETURN';
  170. '      END';]
  171.  
  172. jac22=[...
  173. '      SUBROUTINE JAC22 (T, Y, ydot, PD, CJ, RPAR, IPAR)';
  174. ' ';
  175. '      IMPLICIT DOUBLE PRECISION (A-H,O-Z)';
  176. '      INTEGER  NROWPD';
  177. '      DOUBLE PRECISION T, Y, PD';
  178. '      PARAMETER (NROWPD=2)';
  179. '      DIMENSION Y(2), PD(NROWPD,2)';
  180. '      PD(1,1) = 0.0D0';
  181. '      PD(1,2) = 1.0D0';
  182. '      PD(2,1) = -200.0D0*Y(1)*Y(2) - 1.0D0';
  183. '      PD(2,2) = 100.0D0*(1.0D0 - Y(1)*Y(1))';
  184. '      PD(1,1) = CJ - PD(1,1)';
  185. '      PD(1,2) =    - PD(1,2)';
  186. '      PD(2,1) =    - PD(2,1)';
  187. '      PD(2,2) = CJ - PD(2,2)';
  188. '      RETURN';
  189. '      END';]
  190.  
  191.  
  192. gr22=[...
  193. '      SUBROUTINE GR22 (NEQ, T, Y, NG, GROOT, RPAR, IPAR)';
  194. '      IMPLICIT DOUBLE PRECISION (A-H,O-Z)';
  195. '      INTEGER NEQ, NG';
  196. '      DOUBLE PRECISION T, Y, GROOT';
  197. '      DIMENSION Y(*), GROOT(*)';
  198. '      GROOT(1) = Y(1)';
  199. '      RETURN';
  200. '      END';]
  201.  
  202. //Uncomment lines below: link may be machine dependent if some f77 libs are 
  203. //needed for linking
  204. //unix_g('cd /tmp;rm -f /tmp/res22.f');unix_g('cd /tmp;rm -f /tmp/gr22.f');
  205. //unix_g('cd /tmp;rm -f /tmp/jac22.f');
  206. //write('/tmp/res22.f',res22);write('/tmp/gr22.f',gr22);write('/tmp/jac22.f',jac22)
  207. //unix_g("cd /tmp;make /tmp/res22.o");unix_g('cd /tmp;make /tmp/gr22.o');
  208. //unix_g('cd /tmp;make /tmp/jac22.o');
  209. //          2  Linking the routines
  210. //link('/tmp/res22.o','res22');link('/tmp/jac22.o','jac22');link('/tmp/gr22.o','gr22')
  211. //rtol=[1.d-6;1.d-6];atol=[1.d-6;1.d-4];
  212. //t0=0;y0=[2;0];y0d=[0;-2];t=[20:20:200];ng=1;
  213. //info=list([],0,[],[],[],0,0);
  214. //          3 Calling the routines by dasrt
  215. //[yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,'res22','jac22',ng,'gr22',info);
  216. // hot restart
  217. //[yy,nn,hot]=dasrt([y0,y0d],t0,t,atol,rtol,'res22','jac22',ng,'gr22',info);
  218. //t01=nn(1);t=100:20:200;[pp,qq]=size(yy);y01=yy(2:3,qq);y0d1=yy(3:4,qq);
  219. //[yy,nn,hot]=dasrt([y01,y0d1],t01,t,atol,rtol,'res22','jac22',ng,'gr22',info,hot);
  220.  
  221.  
  222.